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Motivated by the possibility of superconductivity in doped graphene sheets, we investigate superconducting 
order in the extended Hubbard model on the two-dimensional graphene lattice using the variational cluster ap¬ 
proximation (VGA) and the cellular dynamical mean-field theory (CDMFT) with an exact diagonalization solver 
at zero temperature. The nearest-neighbor interaction is treated using a mean-field decoupling between clusters. 
We compare different pairing symmetries, singlet and triplet, based on short-range pairing. VGA simulations 
show that the real (nonchiral), triplet p-wave symmetry is favored for small V, small on-site interaction U or 
large doping, whereas the chiral combination p + /p is favored for larger values of V, stronger on-site interac¬ 
tion U or smaller doping. GDMFT simulations confirm the stability of the p-\-ip solution, even at half-filling. 
Singlet superconductivity (extended 5-wave or J-wave) is either absent or sub-dominant. 


I. INTRODUCTION 

Following the production of graphene sheets in 2004 [1], 
many have speculated about the theoretical possibility of su¬ 
perconductivity in doped graphene. The theoretical discussion 
has been enlarged to include models of interacting electrons 
on the graphene (honeycomb) lattice, without necessarily fo¬ 
cusing on parameter values relevant to graphene, as other sys¬ 
tems based on this geometry exist: for a recent review see, 
e.g., Ref. [2]. A constant source of excitement is the gen¬ 
eral conclusion that superconductivity, if it occurs, should be 
chiral, i.e., should break time-reversal symmetry. This im¬ 
plies the possibility of unidirectional transport along the sam¬ 
ple edge and, with the added effect of spin-orbit coupling, the 
presence of Majorana fermions along the edge, with poten¬ 
tial applications to robust quantum computing. In particular, 
certain vortex excitations mp^ip superconductors have zero- 
energy Majorana modes [3] in their cores, which endow these 
vortices with non-Abelian statistics [4]. 

Many studies predict a superconducting order parameter 
with d + id symmetry, i.e., a chiral state based on singlet pair¬ 
ing [5-10]. But many of those have excluded triplet pairing 
from the outset. A recent Quantum Monte Carlo study [11] 
compares singlet and triplet states and predicts that the favored 
state would have p^ ip (triplet) symmetry [12], but it is per¬ 
formed at low density (^ 20%), whereas we are interested in 
the vicinity of half filling. Such comparisons between triplet 
or singlet pairing have also been made using other methods or 
different theoretical models [13, 14] and the favored pairing 
state is not strikingly obvious. 

In this work, we add the perspective of two different meth¬ 
ods: the variational cluster approximation (VGA) [15-18] and 
the cellular dynamical mean field theory (CDMFT) [19-22]. 
General reviews on dynamical quantum cluster methods can 
also be found in Refs [23-25]. VCA is a dynamical vari¬ 
ational method: it identifies the best possible electron self¬ 
energy T,{co) in a. restricted space of self-energies that are the 
physical self-energies of a small cluster of atoms. CDMFT 
is based on the same principle, except that the best possible 
self-energy is determined self-consistently in an Anderson im¬ 
purity problem where the small cluster plays the role of the 


impurity. We apply these methods to the extended Hubbard 
model on the graphene lattice, from half-filling to about 20% 
doping and compare various pairings: singlet and triplet, chi¬ 
ral and non chiral. The dominant pairing symmetry found by 
VCA, i.e. the one with the smallest free energy at zero tem¬ 
perature, is p 4-ip: a chiral, spin-triplet pairing. This solution 
is also found using CDMFT. The range of parameters stud¬ 
ied, in particular for the on-site repulsion U and the nearest- 
neighbor repulsion V, contains accepted values for graphene 
sheets. Both methods find that the strength of p 4-ip super¬ 
conductivity increases with V. 

The paper is organized as follow: in Sect. II, we define the 
model and the various pairing symmetries from mean-field the 
point of view. In Sect. Ill we review the VCA and its applica¬ 
tion to systems with extended interactions, before presenting 
our results in Sect. IV. In Sect. V, the CDMFT is applied to 
the same problem. A discussion follows in Sect. VI. 

II. MODEL AND MEAN-FIELD DESCRIPTION 
A. The model 

We consider the extended, one-band Hubbard model de¬ 
fined on the graphene (or honeycomb) lattice, which contains 
two sublattices A and B as illustrated in Fig 1 . The Hamilto¬ 
nian can be expressed as 

H = -t Y, + “ mL"*- 

reA,Gj r 

+ Y, '^r^r+ey (1) 

r reAj- 

where destroys (creates) an electron of spin cr in a Wan- 

nier orbital at site r, rir^o = <^r!a^r ,<7 is the number of electrons 
of spin (7 at site r, and rir = ^r,t + The three vectors 
01 ^ 2,3 link a site of sublattice A with its three nearest neigh¬ 
bors (NN) on sublattice B, and are oriented at 120° of each 
other. The first and last sums run over sites of the A sublattice 
only and contain respectively all hopping terms and extended 
interactions. The other sums run over all sites: p is the chem¬ 
ical potential and U the on-site repulsion. 
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FIG. 1. (Color online) Tiling of the graphene lattice by 6-site clusters 
(gray shading) used in VGA. The A and B sublattices are indicated, 
as well as the three elementary vectors ei 2 , 3 . 


This model constitutes an approximate description of 
graphene sheets, wherein longer-range Coulomb interactions 
and hopping are neglected. In pure graphene the band is half- 
filled and t is estimated at 2.8 eV, the on-site Coulomb repul¬ 
sion U is expected to be around 9.3 eV « 3.3r and the nearest- 
neighbor Coulomb repulsion at V = 5.5 eV ^ 2r [26]. In the 
rest of this paper, we will set t = 1, thus defining the unit of 
energy, and we will work at zero temperature. We are con¬ 
cerned with electronic degrees of freedom only, and neglect 
all phonon-related effects. 


B. Pairing symmetries 

We start by providing a description of superconducting 
(SC) order in the mean-field picture. Ref. [13] provides a com¬ 
plete classification of pairing operators in terms of low-degree 
polynomials of momentum, according to the irreducible rep¬ 
resentations of in momentum space. Here we will fol¬ 
low a different approach, based on a real-space description 
of pairing operators defined on adjacent sites; we neglect the 
possibility of on-site (singlet) pairing because of the on-site 
repulsion U. In momentum space, this amounts to using basis 
functions that are aware of the Brillouin zone, i.e., complex 
exponentials of wavevectors times nearest-neighbor vectors. 
The relevant pairing operators defined on the links between 
adjacent sites are 

singlet: Si r = <^r,t^r+e/,f ~ ^r+e/,t 

triplet: Tfr = Cr,tCr+e,,; + Cr,;^^r+e„t 

Given the three elementary directions on the graphene lattice, 
this makes a total of six operators per site, which can be com¬ 
bined into operators of well-defined symmetry as 


^singlet — ^(Ai^i^r + ^2^52,r + ^35'3,r) 

r 

^triplet = ^(AiTi^r + A272,r + A3r3^r) 


( 3 ) 


TABLE I. Irreducible representations (irreps) of associated with 
the six pairing operators defined on nearest-neighbor sites. Sj and Tj 
are the singlet and triplet pairing along the directions Cy indicated on 
Fig. 1. The last four rows show the chiral representations, which are 
complex combinations of the real operators defined under Ei and E 2 . 


Irrep 

symbol 

operators 

Ai 


A^ = ^ (5l,r + 5'2,r + 5'3^r) 

r 

5l 

/ 

^/ = L(^l,r + ^ 2 ,r + 73,r) 

El 

d 

A(i,l = ^ (‘^hr ~ ‘^ 2 ,r) 



A ^/,2 = l^(5'l,r-5'3,r) 

r 

El 

P 

Vl=I(ri,r-r 2 ,r) 

r 

V 2 = E(7’i,r-73,r) 

r 

chiral representations 

El d + id Ad+id = E (5l,r + 


d — id 

K-id = I {Sl,r + 

El 

PEip 

K+ip = E {Ti ,r + Ti^r + 73 ,r) 


p-ip 

K-iP = E 

r V / 


where the relative amplitudes (Ai,A 2 ,A 3 ) define the symme¬ 
try of each operator. These may be classified according to the 
irreducible representations (irreps) of D^k (or C^v, which is 
equivalent for a purely two-dimensional system). These are 
given in Table I, and illustrated on Fig. 2. Note that the irreps 
A 2 and B 2 do not exist in this six-dimensional space of op¬ 
erators. Representations Ei and E 2 are two-dimensional, and 
only one of their components is illustrated on Fig. 2. This al¬ 
lows for the existence of complex representations d ± id and 
pEip expressed in the last four rows of Table 1. 


C. Mean-field description 

The goal of this subsection is to develop a physical or geo¬ 
metric sense for the superconducting order parameter through 
the mean-field description. We are not performing self- 
consistent mean-field computations, which are not possible in 
the framework of the Hubbard model. The more powerful 
variational cluster approximation (VGA) and cellular dynam¬ 
ical mean-field theory (CDMFT) will be used instead, in the 
next sections. 

The BCS Hamiltonian in momentum space is expressed in a 
Nambu description by arranging the creation and annihilation 
operators for the two sublattices and the two spins into a four- 
component object: 


-k,f) 


( 4 ) 
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FIG. 2. (Color online) Representation of singlet and triplet pairing 
amplitudes in real space. Blue means positive, red means negative. 
Only one component out of two is illustrated for the E\ and E 2 rep¬ 
resentations: 2 ^nd 2- Upon rotating by 60 degrees, the two 

sublattices A and B are interchanged, which changes the sign of the 
triplet pair (because of the anticommutation relations) but not that of 
the singlet pair. 


and 

Bk = + + |j?kP|r7-kP - 2Re(^T7kJ?_kT^) 

+ M^(|T]kP + |J7-kp-2|T1cP) ( 11 ) 

If the amplitudes At are real, i.e., for the real representations 
V, p, d and /, then r 7 _k = but this is not true for the chiral 
representations. In the normal case (T 7 k = 0), Eq. (9) reduces 
to the graphene dispersion E\^ = p±\yi^\. 

The Dirac points K = (2;r/3,2;r/3A/3) and K' = 
(2;r/3, —In 13^/3), located on the Brillouin zone corners, are 
special, since /(K) = 7 (K') = 0. At half-filling (p = 0), the 
normal-state dispersion vanishes at these points and the low- 
energy dispersion is made of cones issueing from them. The 
gap function r]k also vanishes at the Dirac points if all three 
amplitudes A/ have the same phase, i.e. for v and / symme¬ 
tries. This implies that v-wave and /-wave superconductiv¬ 
ity is gapless at half filling. So is chiral superconductivity 
(d + id and p-\-ip). Indeed, in the right-handed case ( 8 ) one 
has r 7 K' = 0 and r 7 _K' 7 ^ 0. Nevertheless, ^k' = 0 at half-filling 
and the gap vanishes. The same is true at the other Dirac point, 
since t]k 7 ^ 0 and r]_K = 0. Thus, superconductivity may be 
“hidden”, or gapless at half-filling [10]. 


The Hamiltonian then takes the form: 


k 

with the 4 X 4 Hermitian matrix: 


/-M 0 -r/A 

t Cn-k 0 

0 CCk M 
\-j?k 0 n / 


( 6 ) 


where p is the chemical potential and where 

^ = -/ £ and 77k = £ (7) 

7 = 1 , 2 ,3 7 = 1 , 2 ,3 


The symbol in ( 6 ) is 1 and — 1 for triplet and singlet pairing, 
respectively. The pairing amplitudes At may be read from the 
definition of the operators of Table I. For instance, in the case 
of 2 i p ^ ip or d id symmetry, they are 

(Ai,A2,A3)=A(l,e2'^'/3,e4^'/3) (8) 

where A is a global pairing amplitude. 

It is not difficult to show that the dispersion relation derived 
from the mean-field Hamiltonian ( 6 ) is 

£’k = (9) 

where 

^k = M^ + l)tP + ^|J7kP + ^l»?-kP (10) 


D. Order parameter 

The momentum-dependent superconducting order param¬ 
eter is defined as the integral over frequency of the 

Gorkov function (the anomalous part of the Green function): 

'¥ab{k)= r ^Fab{k,i(0) ( 12 ) 

j —oo 27r 

Here {a^b) are sublattice indices taking two possible values. 
The Gorkov function Fab is the top-right block of the Nambu 
Green function G^v(k, (o) defined as follows at zero tempera¬ 
ture: 


+ (QK(k) ^^^_^^ -4(k)|n) , (13) 

where m is a complex-valued frequency, |^ 2 ) is the many-body 
ground state and £"0 the ground state energy. The indices p , V 
take the four possible values defined in (4). For a two-band 
model, Fab = Ga^b+2- In the special case of the non-interacting 
BCS Hamiltonian (5), the Nambu Green function is 

and the order parameter ^^ab (k) is readily calculated from the 
negative eigenvalues of and the corresponding eigenvec¬ 
tors. We illustrate on the left panel of Fig. 3 the supercon¬ 
ducting order parameter of type p^ ip with A = 0.3 at 10% 
doping. Note how the phase of the order parameter An circles 
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FIG. 3. (Color online) Color representation of the superconducting order parameter as a function of wave vector, for a symmetry. 

Left half of the figure: a BCS state; right half: a solution found in VCA. Panels labeled (a) and (c) represent the amplitude (red is maximum, 
blue means zero). Panels labeled (b) and (d) represent the phase (the phase color map is shown in the middle; white means an amplitude 
lower than some cutoff value). The top panels ((a) and (b)) represent the intra-sublattice component 'Pn(k). The bottom panels ((c) and (d)) 
represent the inter-sublattice component 'Pi 2 (k). The Brillouin zone is indicated, as well as the Fermi surface around one of the Dirac points 
K (dashed curve). Doping was set at 10%. The VCA solution was obtained at (C/, V) = (3,0.4). 


once around the Dirac points K and K'. By contrast, the phase 
of Ai 2 circles twice around K and does not circle around K', 
whereas the opposite would be true for p —ip ox for A 21 . 

It is difficult to encapsulate the “amount” of SC order in 
a simple number. The best choice for that is the root-mean- 
square SC order parameter: 

= (15) 

abJ 

This definition has the advantage of being invariant under 
changes of basis affecting the sublattice (or band) indices. 


III. THE VARIATIONAL CLUSTER APPROXIMATION 
A. The method 

The variational cluster approximation (VCA) [15] can be 
used to investigate the possibility of superconductivity in 
Model (1). VCA has been used extensively, in particular to 
study the emergence of J-wave superconductivity in a sim¬ 
ple description of the high-Tc cuprates based on the square- 
lattice, one-band Hubbard model [16, 17]. It is based on Pot- 
thoff’s self-energy functional approach [27] (for a review, see 
Ref. 28). In VCA, the lattice is tiled into an infinite collection 
of identical clusters, like the six-site cluster used in this work 
and illustrated on Fig. 1. We must distinguish the original 


Hamiltonian H, defined on the infinite lattice, from a reference 
Hamiltonian obtained from H by (1) severing the inter¬ 
cluster hopping terms and (2) adding a small number of Weiss 
field in order to probe certain broken symmetries. Any one- 
body term can also be added to H' \ in particular, the chemical 
potential p' of H' may be different from the one appearing in 
H. The basic requirement is that H and H' share the same in¬ 
teraction term. Even though H' is defined on an infinite set of 
disconnected clusters, in practice we work on its restriction to 
a single cluster, since all clusters are identical. 

The self-energy T,{co) associated with H' is used as a vari¬ 
ational self-energy, in order to construct the Potthoff self¬ 
energy functional: 

+ Trln[-(Go' -S(/j))-'] - Trln(-G'(/z)) (16) 

where G' is the physical Green function of the reference sys¬ 
tem, Go is the non-interacting Green function of the origi¬ 
nal system and h denotes collectively the coefficients of all 
the adjustable one-body terms added to H'. The symbol Tr 
stands for a functional trace, i.e., a sum over frequencies, mo¬ 
menta and bands. Q! is the ground state energy (chemical 
potential included) of the cluster which, along with the as¬ 
sociated Green function G', is computed numerically, in our 
case via the exact diagonalization method at zero temperature. 
Eq. (16) provides us with an exact, nonperturbative value of 
the Potthoff functional ^l[5](/z)], albeit on a restricted space 
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of self-energies '^{h) which are the physical self-energies of 
the reference Hamiltonian H'. Expression (16) is computed 
numerically in order to look for stationary points of that func¬ 
tional, for instance via Newton’s method. The resulting value 
of h defines the best possible self-energy H for that parame¬ 
ter set; it is then combined with Go to form an approximate 
Green function G for the original Hamiltonian H, from which 
any one-body quantity, for instance the order-parameters as¬ 
sociated with broken symmetries, can be computed. 

When confronted with competing solutions, obtained for 
instance via different sets of Weiss fields, the one with the 
lowest value of the Potthoff functional is selected. VGA re¬ 
tains the correlated character of the model, since the local in¬ 
teraction is not factorized in any way. The approximation may 
be controlled in principle by varying the size of the cluster and 
the number of variational parameters used. 

Since one of the goals of this work is to identify the sym¬ 
metry of the superconducting order parameter in Model (1), 
we will use the 6-site, ring cluster depicted on Fig. 1, because 
it is the most symmetric we can use, even though it is not the 
largest. A ten-site cluster will also be used in order to as¬ 
sess the robustness of our predictions. For all calculations in¬ 
volving superconductivity, we used the Nambu formalism, in 
which a particle-hole transformation is performed on the spin- 
down orbitals. The pairing operators then have the appearance 
of hopping terms; two of the three components of the triplet 
pairing operator cannot be easily described that way, but rota¬ 
tion invariance allows us to concentrate on the 5^ = 0 compo¬ 
nent. 


B. Extended interactions 

The VGA approximation as summarized above only applies 
to systems with on-site interactions, since the Hamiltonians H 
and H' must differ only by one-body terms, i.e., they must 
have the same interaction part. This is not true if extended 
interactions are present, as they are truncated when the lattice 
is tiled into clusters. To treat the extended Hubbard model, 
one must, in addition, perform a Hartree approximation on 
the extended interactions. We call this the dynamical Hartree 
approximation (DHA), as the on-site and extended interac¬ 
tions are treated exactly within the cluster. It has been used 
in Ref. 29 in order to assess the effect of extended interactions 
on strongly-correlated superconductivity. 

The extended interaction in the model Hamiltonian (1) is 
replaced by 

\ ^ + nr«r' “ Wr^r') (17) 

r,r' r,r' 

where denotes the extended interaction between sites be¬ 
longing to the same cluster, whereas those interactions 
between sites belonging to different clusters. Both the first 
term (V^) and the second term which is a one-body op¬ 
erator, are part of the lattice Hamiltonian H and of the cluster 
Hamiltonian H'. The mean fields fir must be determined self- 
consistently via a repeated application of the VGA method or 



U 

FIG. 4. (Color online) Half-filling phase diagram of the C/ — V ex¬ 
tended Hubbard model, obtained through VGA. The critical U for 
the appearance of antiferromagnetism at V = 0 is = 3.86. The 
three phases are a charge-density wave (CDW), a normal semi-metal 
(N) and an antiferromagnet (AF). 


even, more simply, of cluster perturbation theory (CPT) [30]. 

In practice, the symmetric matrix is diagonalized and 
the problem is expressed in terms of eigenoperators : 


V" = 


1 






( 18 ) 


For the 6-site cluster used in this work, the two eigenopera¬ 
tors considered correspond to a uniform shift of the chemical 
potential and a charge-density-wave operator: 

»Il=E«r OT2 = E”r-E"*- (1^) 

r reA reB 


with the appropriate coupling constants Di = ^ and D 2 = 


IV. RESULTS FROM THE VGA 
A. Antiferromagnetism and charge order at half-filling 

It is well established that the Hubbard model on the 
graphene lattice with nearest-neighbor hopping has an antifer¬ 
romagnetic ground state beyond a certain value of the on-site 
interaction 1/ at V = 0 [31-34]. We mapped the antiferromag¬ 
netic transition line on the U — V plane at half-filling, with 
VGA, in the dynamical Hartree approximation described in 
Sect. IIIB with the mean field m\ of Eq. (19). The Weiss field 
added to the cluster Hamiltonian in order to probe antiferro¬ 
magnetism in VGA is 


^AF = ^AF \ ^ ~ f • 

[ i^A ieB ) 

The phase boundary between the antiferromagnetic (AF) and 
normal phases found in this way is indicated in green on 
Fig. 4. This is a continuous (second-order) transition. At V = 
0, the critical value of U for antiferromagnetism is Uc ^3.86, 
identical to 2 decimal places to the value obtained from large- 
scale quantum Monte Carlo simulations [32], and close to 
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Uc ^3.6 from the dynamical cluster approximation [35]. The 
antiferromagnetic phase is bound to extend somewhat away 
from half-filling. Ideally we would want to avoid clashing 
with this phase in our study of superconductivity. 

At larger values of V a charge density wave (CDW) insta¬ 
bility is bound to occur. We determined the phase boundary 
between the CDW and normal phases by applying the dynam¬ 
ical Hartree approximation with the mean fields mi and m 2 of 
Eq. (19), and by comparing the energies Q, of the two com¬ 
peting solutions: a normal state (NS) with m 2 = 0 and a CDW 
with m 2 7 ^ 0. No Weiss field was added in this case, in order 
to simplify the computation: thus CPT was used instead of 
VC A, but the functional (16) was computed as a best estimate 
of the free energy. The phase boundary is indicated in blue 
on Fig. 4, and tends asymptotically towards the line V = U/3 
(dashed line), as expected in the strong coupling limit. For 
sufficiently large U, the two phase boundaries (AF-NS and 
CDW-NS) will cross and a competition between CDW and 
AF phases would need to be examined. The NS-CDW phase 
boundary is basically identical with the one found with the dy¬ 
namical cluster approximation (DCA) on large clusters [35]. 
At f/ = 0 the NS-CDW transition is continuous, but becomes 
discontinuous beyond some value of U. This is in fact an im¬ 
portant test of the dynamical Hartree approximation used in 
the rest of this work. In particular, the constant correction 
added to the energy (the last term of Eq. (17)) is crucial if the 
phase boundary is to tend towards the line V = f//3. 

Curiously, the values (U^V) = (3.3,2.0) computed in 
Ref. [26] lie within the CDW phase. But adding a second- 
neighbor Coulomb interaction V' would push the CDW phase 
boundary further up: the phase boundary in the strong¬ 
coupling limit is easily seen to be the line V = t//3 4-2y', 
and the value of V' computed in Ref. [26] for graphene is more 
than enough to push that phase boundary beyond the proposed 
values of {U,V). 


B. Superconductivity 
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FIG. 5. (Color online) Condensation energy in units of t of the var¬ 
ious triplet and singlet superconducting solutions found in VGA as 
a function of doping d = I —n. Top left: Stable solutions at C/ = 3 
and y = 0. Bottom left: Stable solutions at C/ = 3 and V = 0.4. Top 
right: Stable solutions aiU = 6 and y = 0. Bottom right, the same, 
but computed on a larger, 10-site cluster made of two edge-sharing 
hexagons. The J-wave solution was not stable, and the >s-wave solu¬ 
tion was only stable for (C/, y) = (3,0.4) among the solutions found, 
and was never the lowest-energy solution. 
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- V = 0.2 - V = 0.4 - V = 0.6 
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In VGA, the possible presence of superconductivity is 
probed by adding to the cluster Hamiltonian H' one of the 
pairing operators appearing in Table I. The VGA computa¬ 
tions for superconductivity use two Weiss fields: the over¬ 
all pairing amplitude A and the cluster chemical potential jd'. 
Treating the latter as a variational parameter guarantees ther¬ 
modynamic consistency [36]. These computations are carried 
for different pairing symmetries, and, when in the presence of 
extended interactions, by performing an extra self-consistency 
loop for the cluster Hartree approximation, as described in 
Sect. Ill B above. 

In order to compare solutions obtained with Weiss fields of 
different pairing symmetries, we compute their energy densi¬ 
ties E = Q,-\- fin, as a function of electron density n. More 
precisely, we compare their condensation energies F’n — Esc, 
the difference between the energy of the normal solution, ob¬ 
tained by using only the chemical potential /i' as a variational 
parameter, and that of the superconducting solution, obtained 
by using both /i' and A as variational parameters. Figure 5 


FIG. 6. (Color online) Left panel: Condensation energy in units of 
t for the preferred superconducting solutions at C/ = 3 as a function 
of doping d, for different values of V. Dashed curves: p-wave solu¬ 
tions; full curves: pE ip solutions. Right panel: Phase diagram on 
the y — d plane at C/ = 3. There is a transition between p and p-\-ip 
solutions. 


shows the condensation energy as a function of doping for 
(f/,y) = (3,0), (3,0.4) and (6,0). In the latter case, two clus¬ 
ter sizes (6 and 10 sites) were used (see below). Different 
pairing symmetries were studied, but triplet pairing is domi¬ 
nant always. The value U = 3 on the left panels was chosen 
because of the absence of antiferromagnetism. 

The first striking result is the existence of triplet-pairing so¬ 
lutions (/, p and p^ ip) even at half-filling. Singlet-pairing 
solutions do not exist at y = 0 , but arise in the presence of 
the extended interaction. However, their condensation energy 
is either too small to appear on the graph, or is smaller than 
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that of triplet-pairing solutions. In particular, a J-wave so¬ 
lution was not found: the Potthoff functional was stationary 
only for a vanishing value of the corresponding Weiss field. 
A chiral, d + id solution is found at (f/, V) = (3,0.4), but its 
condensation energy is negligible and would be barely visible 
if shown on Fig. 5. An extended 5--wave solution is found at 
(t/,y) = (3,0.4), but is never the most stable solution. That 
title goes to the p-wave or to the chiral ip solution, de¬ 
pending on doping and on V. At (f/,y) = (3,0), the dom¬ 
inant solution has p-wave symmetry, but already at y = 0.4 
the chiral, p^ ip solution dominates. 

The left panel of Fig. 6 shows the condensation energy for 
the lowest-energy solution as a function of doping for different 
values of the extended interaction V, at U = 3. The lowest two 
values of y (0 and 0.2) prefer a real, /7-wave solution, whereas 
higher values of V favor the p-\-ip solution. According to 
these results, V has a favorable effect on superconductivity. 
The right panel of Fig. 6 shows a phase diagram on the y — 5 
plane: lower doping and higher values of V favor the chiral 
p^-ip state compared to the non-chiral /7-wave state. 

We now move to stronger coupling. The two panels on the 
right of Fig. 5 shows the condensation energy for the differ¬ 
ent SC solutions (all triplet) at f/ = 6. In principle the solution 
should be antiferromagnetic for a range of doping around half¬ 
filling for this value of U. Here antiferromagnetism was sup¬ 
pressed in order to simplify computations: we are concerned 
here with the preferred SC pairing, not the possible coexis¬ 
tence with antiferromagnetism. The top right panel shows the 
condensation energy as a function of doping computed from 
the six-site cluster illustrated on Fig. 1 and used in most of 
this work. The bottom panel shows the corresponding results 
on a ten-site cluster made of two hexagonal cells. Using a 
larger cluster provides a check on the robustness of VC A re¬ 
sults, even though a finite-size analysis is rarely possible. Our 
results still stand, except that the 10-site cluster does not have 
the full symmetry of the lattice, and therefore a Weiss field 
of one symmetry (for instance p^ ip) will lead to a nonzero 
average of the /-wave pairing operator as well, which would 
not happen for the hexagonal, 6-site cluster. 


V. CLUSTER DYNAMICAL MEAN-FIELD THEORY 

We use cellular dynamical mean field theory (CDMFT) to 
confirm the appearance of triplet superconductivity by an in¬ 
dependent method. CDMFT, like VC A, proceeds by tiling 
the lattice with clusters and by computing an optimized self¬ 
energy for each cluster. Unlike VC A, the space of self¬ 
energies is not explored by adding Weiss fields on the cluster, 
but rather by coupling each cluster to a bath of uncorrelated, 
auxiliary orbitals that represent the effect of the cluster’s en¬ 
vironment [19-22]. The cluster Hamiltonian is supplemented 
by bath-cluster hybridization and bath energy terms: 

^hyb = ^ +H.C. 

u.a 

t 

^bath — / , ^uv^ii^v 


£3 



FIG. 7. (Color online) Left: Cluster-bath system used in CDMFT. 
The cluster sites are black circles, the bath orbitals red squares. The 
bath parameters are orbital energies £/, hybridizations Ot and intra¬ 
bath triplet pairing di (the arrows indicate the conventional sign of 
that pairing, which is odd under spatial inversion). The chiral super¬ 
current I is measured along the green dotted triangle. The four-site 
cluster by itself does not tile the lattice, but can be combined with its 
inverted image to tile the lattice (right panel). 


where denotes the annihilation operator for the bath orbital 
labeled p. Again, we use the Nambu formalism, wherein a 
particle-hole transformation is applied to the spin-down or¬ 
bitals. Hence the matrices Oai^ and may contain off- 
diagonal blocks associated with pairing, contributing to the 
anomalous Green function. 

The Hamiltonians (21), together with the restriction of the 
Hubbard Hamiltonian (1) to the cluster, defines an Anderson 
impurity model. The cluster Green function, when traced over 
the bath orbitals, takes the following form as a function of 
complex frequency CO\ 

G'-\(0) = (0-t-r((0)-'Eic0) (22) 

where the hybridization matrix T{co) is 

T{co)=e{co-e)-^e'^. (23) 

in terms of the matrices and e^y. In practice, the clus¬ 
ter Green function is computed from an exact diagonalization 
technique using variants of the Lanczos method (just like in 
VGA) and the self-energy is extracted from Eq. (22). 

The Green function G(k, (o) for the lattice model is then 
computed from the cluster’s self-energy as 

G-Hk,fi)) = Go'(k,fi))-S(fi)) (24) 

Here k denotes a reduced wave vector, belonging to the Bril- 
louin zone associated with the superlattice of clusters that de¬ 
fines the tiling, and Go is the non-interacting Green function. 
All Green function-related quantities are 2Nc x 2Nc matrices, 
Nc being the number of sites in the unit cell of the superlattice, 
which is made of one or more distinct clusters (the factor of 2 
is there because of spin, or more precisely Nambu space). 











8 



0 0.02 0.04 0.06 0.08 0.1 0.12 

5 

FIG. 8. (Color online) Top panel: expectation value of the p-\-ip 
pairing operator in the solutions found by CDMFT for C/ = 3 
and different values of the nearest-neighbor Coulomb interaction V, 
as a function of doping. Bottom panel: corresponding value of the 
chiral current I circulating along the triangle defined by the sites 2, 3 
and 4 of the cluster. 

The bath and hybridization parameters {Sj^y^Oap) are de¬ 
termined by the self-consistency condition 

G'(ffl) = y — £G(k,fi)) (25) 

k 

where the integral is carried of the reduced Brillouin zone 
(the domain of k). In other words, the local Green function 
G'(m) should coincide with the Fourier transform of the full 
Green function at the origin of the superlattice. This condi¬ 
tion should hold at all frequencies, which is impossible in a 
zero-temperature implementation of CDMFT because of the 
finite number of bath parameters at our disposal. Therefore, 
condition (25) is only approximately satisfied, through the use 
of a merit function. Details can be found, for instance, in 
Refs [22, 37, 38]. 

In the presence of extended interactions, the dynamical 
Hartree approximation is used in conjunction with CDMFT, 
but in that case the mean-field parameters are converged at the 
same time as the bath parameters, which makes the method 
more efficient than its VGA counterpart. 

In the zero-temperature formalism used here, the size of the 
bath is limited by the use of the Lanczos method to solve the 
Anderson impurity problem. We used the bath-cluster system 
illustrated on Fig. 7. The cluster has four sites. Together with 
its inverted image, it forms a periodically repeated ‘superclus¬ 
ter’ that tiles the lattice. The six bath orbitals are assigned 
fictitious positions that represent the neighboring sites of the 
cluster. They each have a bath energy £i and a hybridization 6 i 
with the cluster. In addition, the superconducting pairing takes 
place only between the bath orbitals, along the links indicated 
in red on Fig. 7. The triplet pairing amplitudes di are defined 
in the order shown, and may be constrained into parameters 
representing various pairing symmetries. For instance, two 


(f) 

independent /-wave bath parameters Jj 2 could be introduced 
such that 

(f) (f) 

Ji = J 3 = J 5 = d\ d2 = d4 = d^= J/ (26) 

whereas two p + //7-wave parameters would be intro¬ 

duced such that 

. ^ (27) 

d2 = e-2'V3^4 = 

and p —ip bath parameters would be defined by complex con¬ 
jugation of the prefactors. Using bath parameters with the 
proper symmetries helps confirming that the converged values 
do indeed represent solutions with well-defined symmetry¬ 
breaking patterns, as some of these parameters, of a given 
symmetry, will converge to nonzero values whereas all oth¬ 
ers will converge to zero. 

In order to facilitate convergence we have set all 6i to a 
common value and arranged the bath energies into two groups: 
£i = £3 = £5 = £ and £2 = £4 = £6 = • In studying p + ip 

pairing, we thus have a total of 3 -b 4 = 7 variational parame¬ 
ters (the 4 come from the real and imaginary parts of 

and plus an optional 6 others if other superconduct¬ 

ing channels are put in competition with the p^ ip channel, 
to check stability. 

Figure 8 shows the results of CDMFT computations per¬ 
formed on the cluster-bath system depicted on Fig. 7 for U = 3 
and several values of the the nearest-neighbor repulsion V, as 
a function of doping 5. The top panel shows the absolute 
value of the expectation value {i^p+ip) of the p^ ip pairing 
operator defined in Table (I). The bottom panel shows the 
chiral current /, defined on Fig. 7, that circulates around the 
cluster, as a function of doping. This is a direct measure of 
the chiral character of superconductivity on the cluster. The 
current actually flows from site 2 to site 3 via the bath sites, 
and so on, in a circular manner, back to site 2. This can only 
be measured on the cluster: physically, it would vanish in the 
bulk and only appear as an edge effect. 

The chiral p + ip solution is indeed found, in preference 
to non-chiral solutions. In other words, when the 6 anoma¬ 
lous bath parameters are allowed to vary in both their real and 
imaginary parts, the p + /p-wave solution is found, the oper¬ 
ator Ap^ip has a nonzero expectation value and the other op¬ 
erators Ap-ip and Af have zero expectation value. Of course, 
initial values of the bath parameters determine whether the 
/7 + //7 or the p —ip solutions emerge in the end. Like in the 
VGA solutions, the nearest-neighbor repulsion V favors su¬ 
perconductivity. 

As the chemical potential jl is varied towards half-filling, 
the p-\-ip order parameter decreases, but does not vanish at 
half-filling. Thus there is a superconducting solution at half 
filling, at the particle-hole symmetric point, like in VGA. 

Similar calculations were attempted for v, d and d id- 
wave superconductivity, with the bath operators dt of Fig. 7 
replaced by singlet pairing operators, but were not successful: 
either the trivial (normal) solution was found, or the CDMFT 
procedure did not converge. 
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VI. DISCUSSION 

Our approach is based on a real-space analysis and is not 
confined to the neighborhood of the Fermi surface. Short- 
range correlations are taken into account exactly, and retarda¬ 
tion effects may be important: we go well beyond mean-field 
theory, even though long-range fluctuations are not taken into 
account. We also study, on the same footing and without bias, 
all possible pairing symmetries. 

Our conclusions are to be contrasted with those of other 
studies that use different methodologies, and that conclude 
that the preferred SC solution has d + id symmetry: For in¬ 
stance, in Ref. [5], a renormalization-group analysis based on 
a small number of k-points and repulsive interactions is per¬ 
formed. In Ref. [7], renormalized mean-field theory is ap¬ 
plied on the r — / model (thus in the strong coupling limit, not 
the same regime as ours). A functional renormalization-group 
analysis is also performed, with the same conclusions, even 
though the approach is usually applied at lower coupling. In 
Ref. [8], a variational Monte Carlo method is applied, but a d- 
wave symmetry is assumed at the outset. The constrained path 
Monte Carlo analysis of Ref. [9] necessitates an initial wave 
function, and thus may be biased towards a singlet-pairing so¬ 
lution as well. 

On the other hand, the mean-field analysis of Ref. [10] con¬ 
cludes that triplet, ip superconductivity is preferred if the 
on-site interaction is repulsive and the nearest-neighbor inter¬ 
action is attractive. But retardation effects from a strong on¬ 
site interaction, not visible in a mean-field treatment, may lead 
to an effective, nearest-neighbor attraction. As demonstrated 
in Ref. [29], an additional, repulsive nearest-neighbor interac¬ 
tion V does not necessarily suppress this effect. This is a com¬ 
plicated dynamical question, especially in the intermediate¬ 
coupling regime. 

Many authors have argued that antiferromagnetic spin fiuc- 
tuations provide the pairing “glue” in high-T^ superconduc¬ 
tors. In particular, within the Hubbard model, the energy scale 


associated with short-range spin fiuctuations has been shown 
to correlate with features of the anomalous self-energy [39]. 
Triplet pairing, on the other hand, would require ferromag¬ 
netic spin fluctuations. 

VII. CONCLUSION 

We have shown that triplet, /7-wave superconductivity 
emerges as the dominant channel for superconductivity in the 
extended Hubbard model on the graphene lattice at weak to 
moderate coupling for dopings ranging from zero to 20%, us¬ 
ing the variational cluster approximation (VGA) and cellular 
dynamical mean field theory (CDMFT) with an exact diago- 
nalization solver at zero temperature. In the presence of an ex¬ 
tended interaction V, we performed a mean-field decoupling 
of the inter-cluster interaction, a method we call the dynam¬ 
ical Hartree approximation (DHA), used in conjunction with 
VGA or GDMFT. 

VGA simulations show that the real (nonchiral) p-wave 
symmetry is favored for small V, small on-site interaction U 
or large doping, whereas the chiral combination p + //7 is fa¬ 
vored for larger values of V, stronger on-site interaction U or 
smaller doping. In this regime, superconductivity exists even 
at half-filling, even though the order weakens on approaching 
half-filling. 

A study of the pairing dynamics similar to that of Refs [29, 
39, 40] would be of interest. 
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